function matpath = read_nd2(nd2_path,matdir,options)

%% Thomas Stevens - Lim Lab
%opened September 2012
%matpath = read_nd2(id,matpath)
%Load .nd2 file, extract metadata, convert to indexable arrays and masks
%Dynamically writes matfile with variables:
%   I = uint16, (x,y,z,t,c,p)
%   BW = logical, (x,y,z,t,c,p)
%   time_v = double, (t,c,p)
%   s = metadata structure
%2D: Uses contrast, 2d weiner, tophat, otsu's method, fill holes, open (despeckle) for masking
%3D: Uses 3D gaussian filter, otsu's method, fill holes, open (despeckle) for masking
%Optional input of custom matfile path, 'default' is:
%~/Desktop/Lim/All_Lim_Microscopy/Yeast/[name].mat
%options = 'Ionly': does not generate binaries
%options = 'par': starts parallel computing on max local cores

% Prompt for a file if not input
if nargin == 0 || exist(nd2_path, 'file') == 0
    [file, path] = uigetfile('*.nd2', 'Select Microscopy File (.nd2)');
    nd2_path = [path file];
    if isequal(path, 0) || isequal(file, 0), return; end
end
if nargin < 2 || (~strcmp(matdir,'default') && exist(matdir, 'dir') == 0)
    matdir = uigetdir('*','Select Save Directory.');
    if isequal(path, 0), return; end
end
if nargin < 3 || isempty(options)
    options = {''};
end

if ischar(options), options = {options}; end

isParallel = startParallel(options);

%% Nested functions for binarization

    function thresh3d(filt_size,top_size,open_size)
        %% 3D Thresholding - can denoise in z as well
        h_bar = waitbar(0,'Filtering/thresholding z-stack timecourse...');
        if isParallel
            close(h_bar)
        end
        G = fspecial3('gaussian',filt_size); %3D Gaussian filter
        H = fspecial3('ellipsoid',top_size); %Ellipsoid structuring element
        if s.z==1
            se = strel('disk',top_size);
        else
            se = strel(H/max(H(:)));
        end
        %Broadcast vars for inner limits of parfor
        [x,y,z,time,channels] = deal(s.x,s.y,s.z,s.time,s.channels);
        channel_name_std = s.channel_name_std;
        otsu_thresh_handle = @otsu_thresh;
        tempDir = [matpath(1:find(matpath=='/',1,'last'))  'workers/'];
        mkdir(tempDir)
        baseFile = [tempDir 'part_%03d.mat'];
        %parfor p=1:s.positions
        parfor p=1:s.positions
            m = matfile(matpath,'Writable',true);
            if isParallel
                %Workers save to distinct files in tempdir
                m_part = matfile(sprintf(baseFile, p),'Writable',true);
            else
                %One core read/write: save directly to mat
                m_part = matfile(matpath,'Writable',true);
            end
            I_p = cell2mat(m.I(p,1));
            BW_cp = false(x,y,z,time);
            for c=1:channels
                disp(['Position #' num2str(p) ', Channel #' num2str(c)])
                if channel_name_std(c)=='W'
                    continue %Skip BF thresholding
                end
                if z==1
                    %Background subtraction -- 2D only (background may be
                    %information in 3D)
                    I_p(:,:,:,:,c) = imtophat(I_p(:,:,:,:,c),se);
                end
                for t=1:time
                    I_mod_t = I_p(:,:,:,t,c);
                    if z==1
                        %Weiner 2d denoising -- 2D only
                        I_mod_t = wiener2(I_mod_t);
                    else
                        %XYZ Gaussian filter -- 3D only
                        I_mod_t = imfilter(I_mod_t,G,'symmetric');
                    end
                    %% Insert 3Ddeconvolution
                    %3D Background subtraction (tophat) -- very slow, no
                    %improvement
                    %m.BW(:,:,:,t,c,p) = imtophat(m.BW(:,:,:,t,c,p),se);
                    %Otsu's method
                    BW_cp(:,:,:,t) = feval(otsu_thresh_handle,I_mod_t,open_size);
                    %Clear all XY borders
                    BW_cp = imclearborder(BW_cp,4);
                    %Clear away transient objects that are not present in
                    %start && end frames
                    if length(time)>1
                        BW_cp = BW_cp &~ imclearborder(BW_cp,time_border_conn);
                    end
                    %Status
                    if ~isParallel && rem(sum([p-1 c-1 t].*[s.time*[s.channels 1] 1]),100)==0
                        waitbar(sum([p-1 c-1 t].*[s.time*[s.channels 1] 1])/prod([s.channels s.positions s.time]),h_bar)
                    end
                end
                m_part.BW(c,p) = {BW_cp};
            end
        end
        if isParallel
            %cat separate worker files
            tic
            disp('Loading temp worker files into central file')
            m = matfile(matpath,'Writable',true);
            for p=1:s.positions
                m_part = matfile(sprintf(baseFile, p));
                for c=1:channels
                    disp(['Position #' num2str(p) ', Channel #' num2str(c)])
                    m.BW(c,p) = m_part.BW(c,p);
                end
            end
            rmdir(tempDir,'s')
            toc
        else
            close(h_bar)
        end
    end

    function I_mod_t = otsu_thresh(I_mod_t,open_size)
        %Otsu's method - pack binary for morphology ops
        thresh = graythresh(I_mod_t);
        I_mod_t = mat2gray(I_mod_t)>thresh;
        %Cancel Flood fill at half full
        if sum(I_mod_t(:)) / numel(I_mod_t) > .5
            I_mod_t = zeros(size(I_mod_t));
        else
            %Uses maximal XYZ connectivity = 26
            %XYZ remove particle if smaller than open_size (cell volume)
            I_mod_t = bwareaopen(I_mod_t,open_size);
            %XYZ hole fill if smaller than open_size (cell volume)
            I_mod_t = ~bwareaopen(~I_mod_t,open_size);
        end
    end

%% Main Function
%Open file
tic
%Default save path for matfile

save_pat = '/([^/]*)\.nd2';
tokens = regexp(nd2_path,save_pat,'tokens');
if strcmp(matdir,'default')
    matpath = ['~/Desktop/Lim/All_Lim_Microscopy/Yeast/' tokens{1}{1} '.mat'];
else
    %user input directory
    matpath = [matdir '/' tokens{1}{1} '.mat'];
end

m = matfile(matpath,'Writable',true);
matpath = m.Properties.Source;
m.initialized = true;
%Only load I if matfile I is empty
if ~ismember('I',who('-file',matpath))
    [label_cell,omeMeta] = bfopen_tvs(nd2_path,matpath);
    s = parse_meta(label_cell,omeMeta);
    m.s = s;
else
    s = m.s;
end
toc

%% Time vector shaping
%Move through time permutations until one is approximately linear --
%metadata is not zct indexed
%if ~ismember('time_v',who('-file',matpath))
tic
disp('Shape time vector')
if s.time==1
    m.time_v = reshape(m.time_v,[s.z s.time s.channels s.positions]);
else
    perm = perms([s.z s.time s.channels s.positions]);
    new_order = 'ZTCP';
    perm_t = perms(new_order);
    std_diff = Inf;
    for p_i=1:length(perm)
        time_test = reshape(m.time_v,perm(p_i,:));
        t_dim = find(perm_t(p_i,:)=='T');
        %lowest mean variation in 1st derivative = most linear
        %should be min second derivative?
        std_diff_test = mean(flat(std(diff(time_test,1,t_dim),[],t_dim)));
        if std_diff_test<std_diff
            std_diff = std_diff_test;
            time_v_linear = time_test;
            old_order = perm_t(p_i,:);
        end
    end
    m.time_v = time_v_linear;
    %rearrange
    for i=1:length(new_order)
        new_dims(i) = find(old_order==new_order(i));
    end
    m.time_v = permute(m.time_v,new_dims);
end
%mean across z-dimension ~ average time through stack
%order = 'ZTCP'->'TCP'
m.time_v = squeeze(mean(m.time_v,1));
toc
%end

%% Initialize and maintain binary as logical, X Y unpacked
m.BW = cell(s.channels,s.positions);
disp('Preallocating false in BW binary...')
tic
for c=1:s.channels
    for p=1:s.positions
        m.BW(c,p) = {false(s.x, s.y, s.z, s.time)};
    end
end
toc

if ~ismember(options,'Ionly')
    %% Pass to nested thresholders
    % um^3 Min Yeast Cytoplasm (to Pixels)
    disp('Threshold')
    tic
    %Transient-clearing connectivity: minimal time only
    time_border_conn = zeros(3,3,3,3);
    time_border_conn(2,2,2,:)=1;
    %Minimum cell size: 20um^3 in voxels
    open_size = ceil(20/prod([s.calibX s.calibY s.calibZ]));
    %1/4 root gaussian filter, calibrated to reflect distance
    filt_size = ceil([s.x^.25 s.y^.25 s.calibX/s.calibZ*s.x^.25]);
    if s.z==1 %2D
        %Volume to area conversion for minimum cell size
        open_size = ceil(open_size^(2/3));
        %sqrt image frame size for background subtraction
        top_size = ceil(sqrt(s.x));
    else
        %3D: calibrated z to scale in distance with x,y
        top_size = ceil([s.x^.5 s.y^.5 s.calibX/s.calibZ*s.x^.5]);
    end
    thresh3d(filt_size,top_size,open_size);
end
toc

end